Read in Analytic Data files
pt.analytic.df <- read.csv(here::here("data", "PupilLightResponse_AnalyticSample_Demog.csv"),
header = T, stringsAsFactors = F)
pt.analytic.df$BMI_c <- pt.analytic.df$BMI - round(mean(pt.analytic.df$BMI), 2)
right_0.post.w <- read.csv(here::here("data", "PupilLightResponse_RightEye_Post_Wide.csv"),
header = T, stringsAsFactors = F)
rownames(right_0.post.w) <- right_0.post.w$subject_id
right_0.post <- read.csv(here::here("data", "PupilLightResponse_RightEye_Post_Long.csv"),
header = T, stringsAsFactors = F)
right_0.post$user_cat2 <- ifelse(right_0.post$user_cat == "daily", 0,
ifelse(right_0.post$user_cat == "occasional", 1, 2))
right_0.post$user_cat2 <- factor(right_0.post$user_cat2,
levels = c(0,1,2),
labels = c("Daily Use", "Occasional Use", "No Use"))
post.user <- ggplot(right_0.post, aes(x=seconds, y = percent_change, group = subject_id))+
geom_line(show.legend = FALSE, alpha = 0.2)+
facet_grid(rows = vars(user_cat2))+
labs(y = "Percent change in pupil response",
x = "Seconds from start of light test")+
theme_bw()
post.user
jpeg(filename = here::here("figs", "Post_use_Indiv_Trajectories.jpeg"),
width = 6, height = 6, res = 300, units = "in")
post.user
dev.off()
Another method that can be used to predict smoking status is a scalar-on-function regression model which uses differences between the whole trajectories to determine smoking status. However, when there are missing values (due to different test durations), all trajectories should be truncated to the shortest or the values should be imputed with fPCA. (See exploration of using SoFR line starting at 3000)
pctChg_vars <- names(right_0.post.w)[grepl("percent_change", names(right_0.post.w))]
pctChg_vars <- pctChg_vars[1:226]
sofr_impute_df <- right_0.post.w[, pctChg_vars]
rownames(sofr_impute_df) <- rownames(right_0.post.w)
sofr_fpca <- fpca.face(as.matrix(sofr_impute_df))
sofr_fpca2 <- sofr_fpca$Yhat
names(sofr_fpca2) <- names(sofr_impute_df)
## Adding in fpca imputed data where missing
sofr_imputed <- sofr_impute_df
sofr_imputed[is.na(sofr_imputed)] <- sofr_fpca2[is.na(sofr_imputed)]
ncols = ncol(sofr_imputed)
pt.analytic.df$female <- ifelse(pt.analytic.df$demo_gender == "Female", 1, 0)
smk.df <- pt.analytic.df[, c("subject_id", "user_cat", "female", "BMI_c")]
smk.df$smoker <- ifelse(pt.analytic.df$user_cat %in% c("daily", "occasional"), 1, 0)
sind <- seq(0, 1, len = ncols)
smat <- matrix(sind, nrow(smk.df), ncols, byrow = T)
# smk.df$Zsm <- I(Zsm)
smk.df$Zraw <- I(as.matrix(sofr_imputed))
smk.df$smat <- I(smat)
smk.df$lmat <- I(matrix(1/ncols, nrow(smk.df), ncols))
smk.df$zlmat <- I(smk.df$lmat * smk.df$Zraw)
# smk.df$bmiMat <- I(replicate(ncol(smk.df$smat), smk.df$BMI))
# smk.df$femaleMat <- I(replicate(ncol(smk.df$smat), smk.df$female))
smk_fglm_ps2 <- gam(smoker ~ s(smat, by=zlmat, bs = "cr", k = 30), data=smk.df,
method = "REML", family = binomial)
summary(smk_fglm_ps2)
##
## Family: binomial
## Link function: logit
##
## Formula:
## smoker ~ s(smat, by = zlmat, bs = "cr", k = 30)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3826 0.5753 2.403 0.0162 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(smat):zlmat 3.195 3.636 11.94 0.0189 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.141 Deviance explained = 13.7%
## -REML = 48.799 Scale est. = 1 n = 84
pt.sofr.roc <- roc(smk_fglm_ps2$model$smoker, smk_fglm_ps2$fitted.values, ci = TRUE)
pt.sofr.roc.auc <- auc(pt.sofr.roc)
pt.sofr.auc.lci <- as.numeric(pt.sofr.roc$ci)[1]
pt.sofr.auc.uci <- as.numeric(pt.sofr.roc$ci)[3]
concurvity(smk_fglm_ps2)
## para s(smat):zlmat
## worst 0.8520874 0.8520874
## observed 0.8520874 0.3376077
## estimate 0.8520874 0.7471882
concurvity(smk_fglm_ps2, full= FALSE)
## $worst
## para s(smat):zlmat
## para 1.0000000 0.8520874
## s(smat):zlmat 0.8520874 1.0000000
##
## $observed
## para s(smat):zlmat
## para 1.0000000 0.3376077
## s(smat):zlmat 0.8520874 1.0000000
##
## $estimate
## para s(smat):zlmat
## para 1.0000000 0.7471882
## s(smat):zlmat 0.8520874 1.0000000
# Sex and gender added as scalar features
smk_fglm_cv <- gam(smoker ~ female + BMI_c +
s(smat, by=zlmat, bs = "cr", k = 30),
data=smk.df,
method = "REML", family = binomial)
summary(smk_fglm_cv)
##
## Family: binomial
## Link function: logit
##
## Formula:
## smoker ~ female + BMI_c + s(smat, by = zlmat, bs = "cr", k = 30)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.00598 0.69602 2.882 0.00395 **
## female -0.87859 0.52508 -1.673 0.09428 .
## BMI_c 0.06205 0.05924 1.047 0.29488
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(smat):zlmat 3.065 3.464 12.02 0.0139 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.168 Deviance explained = 17.2%
## -REML = 48.252 Scale est. = 1 n = 84
anova(smk_fglm_cv, smk_fglm_ps2, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: smoker ~ female + BMI_c + s(smat, by = zlmat, bs = "cr", k = 30)
## Model 2: smoker ~ s(smat, by = zlmat, bs = "cr", k = 30)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 77.138 89.662
## 2 78.922 93.456 -1.7843 -3.7937 0.1251
sigRegion <- function(df, crit.val, imp.var){
time <- imp.var[1]
est <- imp.var[2]
lci <- imp.var[3]
uci <- imp.var[4]
crit.rows <- which((df[, lci] < crit.val & df[, uci] < crit.val) |
(df[, lci] > crit.val & df[, uci] > crit.val))
# print(crit.rows)
group <- vector(mode = "numeric", length = length(crit.rows))
grp.ind <- 1
for(i in 1:length(crit.rows)){
if(i == 1){
group[i] <- grp.ind}
else(if(i > 1){
if(crit.rows[i] - crit.rows[i-1] == 1){
group[i] <- grp.ind
}
else(if(crit.rows[i] - crit.rows[i-1] != 1){
grp.ind <- grp.ind+1
group[i] <- grp.ind
})
})
}
# print(group)
num.grp <- unique(group)
ind.df <- data.frame(group = NA,
start = NA,
end = NA)
for(i in num.grp){
ind.df[i,] <- c(i, min(crit.rows[group == i]), max(crit.rows[group == i]))
}
# print(ind.df)
res.df <- data.frame(int.start = NA,
int.end = NA,
est.time = NA,
est = NA,
lci = NA,
uci = NA)
# print(nrow(ind.df))
for(i in 1:nrow(ind.df)){
sub.df <- df[(ind.df$start[i]):(ind.df$end[i]), ]
# print(sub.df[1, ])
if(sub.df[1, lci] > crit.val & sub.df[1, uci] >crit.val){
start1 <- round(sub.df[1, time], 2)
end1 <- round(sub.df[nrow(sub.df), time], 2)
or1 <- max(sub.df[ , est])
peaktime1 <- round(sub.df[sub.df[, est] == or1, time], 2)
lci1 <- round(sub.df[sub.df[, est] == or1, lci], 2)
uci1 <- round(sub.df[sub.df[, est] == or1, uci], 2)
or1 <- round(or1, 2)
res.df[i, ] <- c(start1, end1, peaktime1, or1, lci1, uci1)
# print(res.df)
}
if(sub.df[1, lci] < crit.val & sub.df[1, uci] < crit.val){
start2 <- round(sub.df[1, time], 2)
end2 <- round(sub.df[nrow(sub.df), time], 2)
or2 <- min(sub.df[, est])
peaktime2 <- round(sub.df[sub.df[, est] == or2, time], 2)
lci2 <- round(sub.df[sub.df[, est] == or2, lci], 2)
uci2 <- round(sub.df[sub.df[, est] == or2, uci], 2)
or2 <- round(or2, 2)
res.df[i, ] <- c(start2, end2, peaktime2, or2, lci2, uci2)
# print(res.df)
}
}
# crit.ind <- by(crit.rows, INDICES = list(group), FUN = function(x) c(min(x), max(x)))
return(res.df)
}
## Covar Adjusted Model
df_sofr_sex <- data.frame("smat" = seq(0, 1, length.out = ncol(smk.df$smat)),
"zlmat" = 1, female = 0, BMI_c = 0)
sofr_pred_sex<- predict(smk_fglm_cv, newdata = df_sofr_sex, se.fit = TRUE, type = "iterms")
# Original Model (unadjusted)
df_sofr <- data.frame("smat" = seq(0, 1, length.out = ncol(smk.df$smat)),
"zlmat" = 1)
sofr_pred<- predict(smk_fglm_ps2, newdata = df_sofr, se.fit = TRUE, type = "iterms")
plot.sofr_sex <- data.frame("frame" = 0:(ncol(smk.df$smat) - 1),
"model_demog" = sofr_pred_sex$fit[,3],
"model_demog_se" = sofr_pred_sex$se.fit[,3],
"model_og" = sofr_pred$fit[,1],
"model_og_se" = sofr_pred$se.fit[,1])
plot.sofr_sex$lci <- exp(plot.sofr_sex$model_demog - 1.96*plot.sofr_sex$model_demog_se)
plot.sofr_sex$uci <- exp(plot.sofr_sex$model_demog + 1.96*plot.sofr_sex$model_demog_se)
plot.sofr_sex$OR <- exp(plot.sofr_sex$model_demog)
plot.sofr_sex$sec <- plot.sofr_sex$frame/30
plot.sofr_sex$p <- exp(plot.sofr_sex$model_demog)/(1+exp(plot.sofr_sex$model_demog))
sofr.sigRegions <- sigRegion(plot.sofr_sex, 1, imp.var = c("sec", "OR", "lci", "uci"))
sofr_plot_sex <- ggplot()+
geom_rect(data = sofr.sigRegions, aes(xmin = int.start, xmax = int.end,
ymin = -Inf, ymax = Inf), fill = "lightpink", alpha = 0.3)+
geom_hline(yintercept = 1, col = "darkorange")+
geom_line(data = plot.sofr_sex, aes(x=sec, y = exp(model_og)))+
geom_line(data = plot.sofr_sex, aes(x=sec, y = exp(model_og-1.96*model_og_se)), linetype = 'longdash')+
geom_line(data = plot.sofr_sex, aes(x=sec, y = exp(model_og+1.96*model_og_se)), linetype = 'longdash')+
geom_line(data = plot.sofr_sex, aes(x = sec, y = exp(model_demog)), col = 'purple')+
geom_line(data = plot.sofr_sex, aes(x=sec, y = exp(model_demog-1.96*model_demog_se)), linetype = 'longdash', col = "purple")+
geom_line(data = plot.sofr_sex, aes(x=sec, y = exp(model_demog+1.96*model_demog_se)), linetype = 'longdash', col = "purple")+
scale_x_continuous(expand = c(0, 0), limits = c(0,7.5),
breaks = c(0,2,4,5,6,7))+
scale_y_continuous(limits = c(0,5),
breaks = c(0,1,2,3,4,5)) +
labs(title = "",
x = "Seconds after start of light test",
y = "Odds ratio between smokers and non-smokers")+
theme_bw()
sofr_plot_sex
right_0.gam.df <- merge(right_0.post, pt.analytic.df,
by = c("subject_id", "user_cat"))
right_0.gam.df$sind <- right_0.gam.df$frame/400
m.post_gam <- mgcv::gam(percent_change ~ s(sind, k=30, bs="cr") +
s(sind, by=occasional, k=30, bs = "cr") +
s(sind, by=daily, k=30, bs = "cr"),
data = right_0.gam.df, method = "REML")
## Create a matrix of residuals
post_gam.resid <- cbind(right_0.gam.df[!(is.na(right_0.gam.df$percent_change)),
c("subject_id", "frame")],
m.post_gam$residuals)
names(post_gam.resid) <- c("subject_id", "frame", "resid")
resid_mat <- reshape(post_gam.resid[, c("subject_id", "frame", "resid")],
timevar = "frame",
idvar = "subject_id",
direction = "wide")
rownames(resid_mat) <- resid_mat$subject_id
resid_mat$subject_id <- NULL
resid_mat <- as.matrix(resid_mat)
post_gam.fpca <- fpca.face(resid_mat, knots = 15)
#Create data frame of eigen functions and attach to original data
Phi_mat <- as.data.frame(post_gam.fpca$efunctions)
colnames(Phi_mat) <- paste0("Phi", seq(1, post_gam.fpca$npc))
Phi_mat$frame <- as.numeric(rownames(Phi_mat)) -1
right_0.gam.df <- merge(right_0.gam.df, Phi_mat,
by = "frame",
all.x = T)
right_0.gam.df <- right_0.gam.df[order(right_0.gam.df$subject_id, right_0.gam.df$frame), ]
right_0.gam.df$subject_id <- as.factor(right_0.gam.df$subject_id)
post_gam.fri <- mgcv::bam(percent_change ~
s(sind, k=30, bs="cr") +
s(sind, by=occasional, k=30, bs = "cr") +
s(sind, by=daily, k=30, bs = "cr") +
s(subject_id, by = Phi1, bs="re") + s(subject_id, by = Phi2, bs="re")+
s(subject_id, by = Phi3, bs="re") + s(subject_id, by = Phi4, bs="re"),
method = "fREML", data = right_0.gam.df, discrete = TRUE)
summary(post_gam.fri)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## percent_change ~ s(sind, k = 30, bs = "cr") + s(sind, by = occasional,
## k = 30, bs = "cr") + s(sind, by = daily, k = 30, bs = "cr") +
## s(subject_id, by = Phi1, bs = "re") + s(subject_id, by = Phi2,
## bs = "re") + s(subject_id, by = Phi3, bs = "re") + s(subject_id,
## by = Phi4, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.3125 0.9197 -11.21 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(sind) 27.54 28.43 115.67 <2e-16 ***
## s(sind):occasional 18.63 21.95 70.14 <2e-16 ***
## s(sind):daily 18.62 21.95 35.42 <2e-16 ***
## s(subject_id):Phi1 82.23 84.00 23472.92 <2e-16 ***
## s(subject_id):Phi2 81.68 84.00 2775.02 <2e-16 ***
## s(subject_id):Phi3 82.23 84.00 887.87 <2e-16 ***
## s(subject_id):Phi4 81.05 84.00 1139.47 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.966 Deviance explained = 96.6%
## fREML = 63334 Scale est. = 2.6314 n = 32642
#original model: post_gam.fri
fosr_sex_gen.fri <- mgcv::bam(percent_change ~
female + BMI_c +
s(sind, k=30, bs="cr") +
s(sind, by=occasional, k=30, bs = "cr") +
s(sind, by=daily, k=30, bs = "cr") +
s(subject_id, by = Phi1, bs="re") + s(subject_id, by = Phi2, bs="re")+
s(subject_id, by = Phi3, bs="re") + s(subject_id, by = Phi4, bs="re"),
method = "fREML", data = right_0.gam.df, discrete = TRUE)
summary(fosr_sex_gen.fri)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## percent_change ~ female + BMI_c + s(sind, k = 30, bs = "cr") +
## s(sind, by = occasional, k = 30, bs = "cr") + s(sind, by = daily,
## k = 30, bs = "cr") + s(subject_id, by = Phi1, bs = "re") +
## s(subject_id, by = Phi2, bs = "re") + s(subject_id, by = Phi3,
## bs = "re") + s(subject_id, by = Phi4, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.117542 0.922950 -10.962 <2e-16 ***
## female -0.206220 0.084028 -2.454 0.0141 *
## BMI_c -0.011684 0.009569 -1.221 0.2221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(sind) 27.54 28.43 115.65 <2e-16 ***
## s(sind):occasional 18.63 21.95 70.40 <2e-16 ***
## s(sind):daily 18.62 21.95 35.37 <2e-16 ***
## s(subject_id):Phi1 82.23 84.00 23642.01 <2e-16 ***
## s(subject_id):Phi2 81.68 84.00 2758.50 <2e-16 ***
## s(subject_id):Phi3 82.23 84.00 884.39 <2e-16 ***
## s(subject_id):Phi4 81.05 84.00 1120.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.966 Deviance explained = 96.6%
## fREML = 63336 Scale est. = 2.6309 n = 32642
anova(fosr_sex_gen.fri, post_gam.fri, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: percent_change ~ female + BMI_c + s(sind, k = 30, bs = "cr") +
## s(sind, by = occasional, k = 30, bs = "cr") + s(sind, by = daily,
## k = 30, bs = "cr") + s(subject_id, by = Phi1, bs = "re") +
## s(subject_id, by = Phi2, bs = "re") + s(subject_id, by = Phi3,
## bs = "re") + s(subject_id, by = Phi4, bs = "re")
## Model 2: percent_change ~ s(sind, k = 30, bs = "cr") + s(sind, by = occasional,
## k = 30, bs = "cr") + s(sind, by = daily, k = 30, bs = "cr") +
## s(subject_id, by = Phi1, bs = "re") + s(subject_id, by = Phi2,
## bs = "re") + s(subject_id, by = Phi3, bs = "re") + s(subject_id,
## by = Phi4, bs = "re")
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32238 84839
## 2 32240 84859 -2.003 -20.631 0.01989 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df_pred_non_d <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 0, "daily" = 0,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0, "female" = 0, "BMI_c" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
df_pred_non <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 0, "daily" = 0,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
pred_f1 <- predict(fosr_sex_gen.fri, newdata=df_pred_non_d, se.fit=TRUE, type = "response")
pred_f2 <- predict(post_gam.fri, newdata=df_pred_non, se.fit=TRUE, type = "response")
df_pred_occ_d <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 1, "daily" = 0,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0, "female" = 0, "BMI_c" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
df_pred_occ <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 1, "daily" = 0,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
pred_f1_1 <- predict(fosr_sex_gen.fri, newdata=df_pred_occ_d, se.fit=TRUE, type = "response")
pred_f2_1 <- predict(post_gam.fri, newdata=df_pred_occ, se.fit=TRUE, type = "response")
df_pred_dly_d <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 0, "daily" = 1,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0, "female" = 0, "BMI_c" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
df_pred_dly <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 0, "daily" = 1,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
pred_f1_2 <- predict(fosr_sex_gen.fri, newdata=df_pred_dly_d, se.fit=TRUE, type = "response")
pred_f2_2 <- predict(post_gam.fri, newdata=df_pred_dly, se.fit=TRUE, type = "response")
df <- data.frame(sec = 0:400/30,
user_group = c(rep("non-user", 802),
rep("occasional", 802),
rep("daily", 802)),
includ_demog = c(rep(1, 401), rep(0, 401),
rep(1, 401), rep(0, 401),
rep(1, 401), rep(0, 401)),
pct_chg = c(pred_f1$fit, pred_f2$fit, pred_f1_1$fit, pred_f2_1$fit, pred_f1_2$fit, pred_f2_2$fit))
df$includ_demog <- factor(df$includ_demog,
levels = c(0,1),
labels = c("Original Model",
"Model - Incl Demographics"))
demog_profiles <- ggplot(data = df, aes(x = sec, y = pct_chg, group = paste0(user_group, "_", includ_demog), col = user_group))+
geom_line(aes(linetype = as.factor(includ_demog)))+
labs(x = "Seconds from start of light test",
y = "Percent change in pupil response",
linetype = "Model Type",
col = "Use Group")
Coefficient plots of demographic model vs original model
df_pred_occ <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 1, "daily" = 0,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0, "female" = 0, "BMI_c" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
df_pred_dly <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 0, "daily" = 1,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0, "female" = 0, "BMI_c" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
pred_f1 <- predict(fosr_sex_gen.fri, newdata=df_pred_occ, se.fit=TRUE, type = "terms")
pred_f2 <- predict(fosr_sex_gen.fri, newdata=df_pred_dly, se.fit=TRUE, type = "terms")
pred_sex_coef_df <- data.frame(seconds = seq(0, 400)/30,
occ_hat = pred_f1$fit[, 4],
occ_se = pred_f1$se.fit[, 4],
dly_hat = pred_f2$fit[, 5],
dly_se = pred_f2$se.fit[, 5])
pred_sex_coef_df$occ_lci <- pred_sex_coef_df$occ_hat - 2*pred_sex_coef_df$occ_se
pred_sex_coef_df$occ_uci <- pred_sex_coef_df$occ_hat + 2*pred_sex_coef_df$occ_se
pred_sex_coef_df$dly_lci <- pred_sex_coef_df$dly_hat - 2*pred_sex_coef_df$dly_se
pred_sex_coef_df$dly_uci <- pred_sex_coef_df$dly_hat + 2*pred_sex_coef_df$dly_se
occ_non.sigRegion <- sigRegion(pred_sex_coef_df, crit.val = 0, imp.var = c("seconds", "occ_hat", "occ_lci", "occ_uci"))
dly_non.sigRegion <- sigRegion(pred_sex_coef_df, crit.val = 0, imp.var = c("seconds", "dly_hat", "dly_lci", "dly_uci"))
post_occ_plt <- ggplot(data = pred_sex_coef_df, aes(x=seconds, y = occ_hat))+
geom_line()+
geom_line(aes(x = seconds, y = occ_hat + 2*occ_se), linetype = "longdash")+
geom_line(aes(x = seconds, y = occ_hat - 2*occ_se), linetype = "longdash")+
# geom_line(aes(x = seconds, y = 0), col = "darkred")+
geom_segment(data=occ_non.sigRegion, aes(x = int.start, y=0,
xend= int.end, yend = 0),
color = "darkred", linewidth =1.2)+
labs(title = "Occasional use vs No use",
y= "",
x = "")+
theme_bw()+
scale_x_continuous(expand = c(0, 0), limits = c(0,10),
breaks = c(0,2,4,6, 8, 10))+
scale_y_continuous(limits = c(-7, 7), breaks = c(-6,-4,-2,0,2,4,6))+
theme(rect = element_rect(fill = "transparent"))
post_dly_plt <- ggplot(data = pred_sex_coef_df, aes(x=seconds, y = dly_hat))+
geom_line()+
geom_line(aes(x = seconds, y = dly_hat + 2*dly_se), linetype = "longdash")+
geom_line(aes(x = seconds, y = dly_hat - 2*dly_se), linetype = "longdash")+
# geom_line(aes(x = seconds, y = 0), col = "darkred")+
geom_segment(data=dly_non.sigRegion, aes(x = int.start, y=0,
xend= int.end, yend = 0),
color = "darkred", linewidth =1.2)+
labs(title = "Daily use vs No use",
y = "",
x = "")+
theme_bw()+
scale_x_continuous(expand = c(0, 0), limits = c(0,10),
breaks = c(0,2,4,6, 8, 10))+
scale_y_continuous(limits = c(-7, 7), breaks = c(-6,-4,-2,0,2,4,6))+
theme(rect = element_rect(fill = "transparent"))
ylab <- textGrob(label = "Difference in percent change in pupil response", rot = 90,
gp = gpar(fontsize = 12))
posval <- textGrob(label = sapply(strwrap("Less Constriction in User", width = 15, simplify = F),
paste, collapse = "\n"), rot = 0,
gp = gpar(fontsize = 8))
negval <- textGrob(label = sapply(strwrap("More Constriction in User", width = 15, simplify = F),
paste, collapse = "\n"), rot = 0,
gp = gpar(fontsize = 8))
post_occ_coef <- grid.arrange(ylab, posval, negval, post_occ_plt,
ncol = 2, nrow = 2,
layout_matrix = rbind(c(1, 2, rep(4, 15)), c(1, 3, rep(4, 15))))
post_dly_coef <- grid.arrange(ylab, posval, negval, post_dly_plt,
ncol = 2, nrow = 2,
layout_matrix = rbind(c(1, 2, rep(4, 15)), c(1, 3, rep(4, 15))))
demog.fosr <- grid.arrange(demog_profiles, post_occ_coef,
post_dly_coef,
layout_matrix = rbind(c(1,2),
c(3,4)))
demog.fosr
## TableGrob (2 x 2) "arrange": 3 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[arrange]
## 3 3 (2-2,1-1) arrange gtable[arrange]